function [y,epsi]=createarp(ar,sigma,nd,nt,cot,pct); 
% function [y,epsi]=createarp(ar,sigma,nd,nt,cot,pct); 
% ar    (1xnlag) autoregressive coefficients
% sigma standard deviation 
% nd    number of obs to discard 
% nt    number of obs to retain 
% cot   if 0 or not specified constant not included 
%       if = 1 constant only 
%       if = 2 constant and trend 
% pct mean or mean and constant coefficient 

randn('seed',sum(100*clock) ); 

ch=rootpol([1 ar]); 
if any( abs(ch) <0.9999 ) 
    warning('Non-stationary model'); 
end 

np=length(ar); 
ng=nd+nt; 
epsi=sigma*randn(ng,1); 

if nargin < 4 
    cot = zeros(ng,1); 
    pct = 0; 
end  

if nargin < 5 | cot == 0; 
    
    cot=zeros(ng,1);
    pct=0; 
    disp('Process without mean'); 
elseif cot==1; 
    
    if length(pct) ~= cot 
        error('Need specify coefficient for constant') 
    end 
    
    cot=ones(ng,1);
    disp('Mean included'); 
elseif cot==2; 
    
    if length(pct) ~= cot 
        error('Need specify coefficient for C and T') 
    end 

    cot=[ones(ng,1) (1:ng)']; 
    disp('Mean and trend included');
end 

pct=r2c(pct); 

y=zeros(ng,1); 
y(1:np)=epsi(1:np); 

ar=ar(:); 
ar_flip = flipud( ar ); 
ar_flip = ar_flip'; 

ii=1; 
for ii=np+1:ng; 
    
    
    y(ii) = cot(ii,:)*pct + ar_flip*y(ii-np:ii-1)  + epsi(ii); 
end; 

y=y(nd+1:end); 
epsi=epsi(nd+1:end); 




